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ABSTRACT 

In this paper we perform a systematic analysis of the structure induced by the onset of gravi- 
tational instabilities in cooling gaseous accretion discs. It is well known that for low enough 
cooling rates the disc reaches a quasi-steady configuration, with the instability saturating at a 
finite amplitude such that the disc is kept close to marginal stability. We analyse the depen- 
dence of the saturation amplitude on the imposed cooling rate, and we find that it scales with 
the inverse square root of the cooling parameter y6 = fcooi/fdyn- This indicates that the heating 
rate induced by the instability is proportional to the energy density of the density waves ex- 
cited by the disc self-gravity. In particular, we find that at saturation the energy dissipated per 
dynamical time by weak shocks due to the gravitational perturbations is of the order of 20 per 
cent of the wave energy. We further perform a Fourier analysis of the disc structure, and sub- 
sequently determine the dominant radial and azimuthal wavenumbers of the density waves. 
While the number of spiral arms (corresponding to the azimuthal wavenumber) is fairly con- 
stant with radius, we find that the disc displays a range of radial wavenumbers whose mean 
increases with increasing radius. The dominant modes closely match the locally most unstable 
wavelength as predicted by linear perturbation analysis. As a consequence, we demonstrate 
numerically that the density waves excited in relatively low mass discs Ma^c/M, ~ 0.1 are 
always close to co-rotation, deviating from it by approximately 10 per cent. This result can 
be understood in terms of the constancy of the Doppler-shifted phase Mach number of the 
flow - the pattern speed self-adjusts so that the flow into spiral arms is always sonic. This has 
profound effects on the degree to which the extraction of energy and angular momentum from 
the mean flow through density waves can be modelled as a viscous process. Our results thus 
provide (a) a detailed description of how the self -regulation mechanism is established for low 
cooling rates, (b) a clarification of the conditions required for describing the transport induced 
by self-gravity through an effective viscosity, (c) an estimate of the maximum amplitude of 
the density perturbation before fragmentation takes place, and finally (d) a simple recipe to 
estimate the density perturbation in different thermal regimes. 

Key words: accretion, accretion discs - galaxies; active - gravitation - hydrodynamics - 
instabilities - planetary systems: formation 



1 INTRODUCTION 

In cold, relatively massive accretion discs the effects of the disc 
self-gravity can become dynamically important, and gravitational 
instabilities may play a major role in determining the long-term 
evolution of the disc. On the one hand the instability can in- 
duce fragmentation of the disc, potentially resu lting in the for- 
matio n of massive stars in AGN accretion discs jNavakshin et al.l 
l2007t) or, in the case of proto -planetary discs, low-m ass compan- 
lons jStamatellosetalJl2007l) or even giant planets ( lBosj|l99l 
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1 19981) . On the other hand, the instability is very effective in 
transporting and redis t ributi ng angular momentum within the disc 
jLodato & Ric3l2004 12005) and might therefore promote accre- 
tion in a variety of diff erent contexts and scales, from supermas- 
sive b lack hole growth ( ShIosman et alj|l990l: Lodato & Natarai^] 
I2OO6I) to star formation jVorobvov & Basull200^ 



The balance between the heating provided through such grav- 
itational instabilities and the cooling as energy is radiated away is 
known to have major implications for the fate of the disc. Frag- 
mentation into bound objects occurs when the cooling dominates, 
whereas the disc can persist in a quasi-equilibrium marginally sta- 
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ble condition if the cooling rate is suffi ciently low ( lGammidl200ll : 
iRice et al.ll20o3 : lLodato & Ricell2004) . 

Such a quasi-steady state is characterised by the presence of 
spiral density waves, which transport energy and angular momen- 
tum through the disc. Where the discs are weakly ionised, and 
therefore where the magneto-rotational instability is likely to be in- 
effective, this angular momentum transport mechanism may be the 
primary driver of accretion. Additionally, these waves extract rota- 
tional energy from the disc, some of which is then returned as heat 
as the waves steepen into shocks. This heat then stabilises the disc 
against further gravitational collapse into bound fragments. As the 
amplitude of the waves increases, so too does the wave energy den- 
sity, increasing the reservoir of energy available t o be returned to 
the disc as heat . However, numerical experiments ( iGaimru 3 1200 ll : 
iRice et al.l2005l) have shown that if the cooling time ?cooi is less than 
a few times the local dynamical time fl"', this feedback process 
breaks down and fragmentation ensues. The quasi-steady marginal 
stability state therefore represents a restricted regime of dynamic 
thermal equilibrium, where the heating through gravitational insta- 
bilities is balanced by the cooling rate. 

In this paper we seek to characterise the relationship between 
the strength of the cooling and the amplitude of the spiral density 
waves excited within the disc through self-gravity, while remain- 
ing within this dynamic thermal equilibrium state. To this end we 
use a Smoothed Particle Hydrodynamics (SPH) code to run global 
numerical simulations of self-gravitating gaseous discs. From such 
controlled numerical experiments, we measure the amplitude of the 
density perturbations over a range of cooling times, down to the 
limit where the disc fragments, and to investigate the overall struc- 
ture formation. Fourier analysis allows us to characterise the mode 
spectra and pattern speeds associated with this structure, and to as- 
sociate the dynamics of the spiral density waves excited through 
self-gravity with the thermodynamics of the disc self-regulation 
process. We begin by discussing some relevant properties of den- 
sity waves and transport processes in gaseous discs in Section 2, 
and link this to the disc thermodynamics in Section 3. Details of 
the numerical simulations and initial conditions are given in Sec- 
tion 4, and the results of these simulations are presented in Section 
5. In Section 6 we discuss the implications of these results and the 
conclusions that may be drawn from them. 



2 DYNAMICS OF SELF-GRAVITATING GASEOUS 
DISCS 

In the following two sections we derive analytically the basic re- 
lations that link the disc quantities in the presence of perturba- 
tions to the induced transport properties of self-gravitating discs. 
These standard results have been mostly derived in the context 
of ste llar dynamics (Shu 1970; Bertin 2000; Binnev & Tremaine 
l2008h . In the context of gaseous (accretion) discs, qualitatively 
similar (but rather more invo lve d) analyses hav e been presented in 
iBalbus & Papaloizoj jl999t) and lBalbuj ( l2003h . 



2.1 Dispersion Relations for Self-Gravitating Discs 

The linear development of the gravitational instability is usually 
described by the dispersion relation D((jj, k, in) between the wave 
(angular) frequency oj and the radial and azimuthal wavenumbers 
of excited modes, k and in respectively. For an infinitesimally thin 
(i.e. 2-dimensional) disc, this can be obtained using the standard 
WKB method of perturbation analysis, and in the tight-winding 



limit (where the radial wavelength is small in comparison to the 
azimuthal wavelength, \kR\ » m, with R the cylindrical radius) we 
have 



D(w, k, m) = (a) - mO.)' 



■ clk^ + 2nGi:\k\ ■ 



V = o, 



(1) 



where C1(R) is the angular speed, and k(R) is the epicyclic fre- 
quency, Cs(R) is the so und speed and £(j?) is the surface density of 
the umperturbed disc jBinnev & Tremainell2008l) . Introducing the 
pattern speed of the wave Q.p(R), where co = mSip, we may express 
this as follows; 



■ Q.f = c^k^ - 2nGI.\k\ + . 



(2) 



Th e well known stability criterion for axisymmetric distur- 
bances ( lToomrelll964) 



Q = ^— > 1 



(3) 



identifies the region of parameter space where the RHS of Eq.Q is 
positive definite, and therefore where the disc is stable at all wave- 
lengths. For the case of an unstable disc, the most unstable wave- 
length fcu„s (i.e. where the RHS of Eq.lO is at a minimum) occurs 
at a radial wave-number 



nOI, 



(4) 



For a marginally stable disc where 2 « I we would therefore ex- 
pect only modes with k » A;;,,,, to be excited, and in this instance, 
Eq.l|2j tells that Qp » SI, i.e. all excited modes are expected to be 
close to co-rotation. Note that the most unstable wave-number is 
exactly the inverse of the disc thickness for a self-gravitating disc 



ttGI. n 



(5) 



where the last approximation is exact for a Keplerian rotation curve. 



2.2 The Stress Tensor 

In non-axisymmetric self-gravitating discs, torques arising due to 
the perturbed gravitational potential play an important role in trans- 
porting both energy and angular momentum. 

Traditio nally, viscous accretion disc s are described by the a- 
formalism of Shakur a & SunvaevI ( Il973l) . where (for an infinitesi- 
mally thin disc) the only non- vanishing component of the vertically 
integrated stress tensor T is the shear term, defined as 



, din n 



(6) 



' dlnR' 

where cf < 1 is a dimensionless parameter which measures the vis- 
cosity. In this case the disc stress is linked to the local pressure Ec^ , 
and indicates that the a-formalism is fundamentally a local rela- 
tionship. Furthermore, since a is not necessarily constant, this rep- 
resents a completely general description for any purely local pro- 
cess. Also note that for Keplerian rotation, dlnQ/dlni? = -3/2, 
implying that the stress is negative, i.e. it acts to oppose rotation, 
and therefore allows for inward accretion flows. 

Accretion disc theory identifies the origin of the stress with 
torques arising due to perturbations in a "turbulent" disc. Crucially 
for the gas dynamics, these perturbations manifest themselves as 
fluctuations in the mean flow velocity, the gravitational potential 
and the magnetic field threading the disc. For non-magnetised self- 
gravitating discs such as we consider here, the stress tensor can 
therefore be broken down into a Reynolds stress term, associated 
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with velocity fluctuations, and a gravitational stress term, associ- 
ated with fluctuations in the gravitational potential. The Reynolds 
stress term r?7" is such that 



(15) 



(V) 



(where Svi{,5v^ are the velocity fluctuations about the mean flow 
velocity in the R and ^ directions respectively and the brackets in- 
dicate azir nuthal averaging) and the gravi tational stress term 
is given bv lLvnden-Bell & KalnaisI ( Il972r) as 



/(IS) 



dz 



(8) 



where again the gR,g^ are the accelerations due to the perturbed 
gravitational potential of the disc in the R and directions respec- 
tively. 

The viscous torque per unit area L„ is related to the vertically 
integrated stress tensor T through 

t„ =RWT, (9) 
which in turn yields 

^" = i^^^^^*) (10) 

oR 

as the only non-zero component of the torque. The power per unit 
surface 5^ produced by this visc ous torque is then given simply by 
jprank et al.ll2002l:|Pringljll98lh 



(11) 



where the subscript a indicates that this relation is expected for a 
viscous disc. Equation (TTJ therefore links the transport of angular 
momentum and the associated rate of work done by torques in the 
case of a local process, as historically modelled via the a-viscosity 
parameter. 

2.3 Wave Energy and Angular Momentum Densities 

In this section we turn our attention to the transport of energy and 
angular momentum through the propagation of the spiral density 
waves expected to arise in a self-gravitating disc. To this end it 
is convenient to introduce the wave action density ^ for density 
waves. Within the WKB approximation used to derive the disper- 
sion relat ion given in Eq.Q, the wave action per u nit surface is 
given by i lToomrelll969l : IShij|l970l ; lFan & Loij|l999h 



m(Q.„ - Q) , 



where (50 is the perturbed gravitational potential, itself given by 
2nG5L 



(12) 



6^ ■■ 



\k\ 



(13) 



where SL is the surface density perturbation. The wave energy per 
unit surface fi„ and the wave angular momentum per unit surface 
-£w are obtained in a straightforward way fro m the wave action 
through the standard wave dynamics relations teertinll200(]| : [shi3 



(14) 



' Note that this system of equations is analogous to that found in quantum 
mechanics for an harmonic oscillator; from the quantum of action fi, the 
quantized energy E and angular momentum S are found via E = ha) and 
5 = tim, where to and m are the angular frequency and spin quantum number 
respectively. 



Combining the first of these relations with Eas.(ll2t and il3l 
we obtain 



SvpVp / SI. 



2 I E / ' 



where 



Vp = inQ.p/k 

Vp = m(Q.„ - Q.)/k 



(16) 



(17) 
(18) 



are the radial phase speed and Doppler- shifted radial phase speed of 
the wave respectively. Note that Eqs.(ll6b and ilSi together explain 
why self-induced density waves are launched at co-rotation - since 
the energy density changes sign at co-rotation, waves that propagate 
away from co-rotation extract no net energy from the flow. 

Looking again at Eas.lll4t and l llSt . we see that the relation- 
ship between the energy per unit surface and the angular momen- 
tum per unit surface in a density wave is given by 



(19) 



In the case of quasi-stationary waves (where Clp is constant) 
propagating in a disc in dynamic thermal equilibrium, the rate at 
which energy is lost per unit surface due to cooling must be bal- 
anced by the power per unit surface £„ dissipated by the waves. 
In order to maintain the amplitude of the wave, the instability has 
to keep extracting energy and angular momentum from the back- 
ground flow. The fluxes of energy (angular momentum) carried by 
the wave are simply £„ (.£«) times the local group velocity. Hence 
when a wave dissipates it adds energy and angular momentum to 
the flow in the ratio of £„ to Jl„, i.e. in the ratio Qp. We therefore 
conclude that 

£„ = fip-Cw. (20) 

Equation J20| ( is the wave analogue of Eq.JTTJ. Comparing 
these two equations, we note a fundamental difference with respect 
to the viscous model - for a given torque waves extract energy 
from the flow at a rate proportional to the wave pattern speed Q.p, 
whereas the rotation speed Q is the underlying rate in the local (vis- 
cous ) case. 

iBalbus & Papaloizoul h999h have noted similarly that in gen- 
eral, energy transport through the gravitational instability cannot be 
described purely in viscous terms, and indeed this is only possible 
at co-rotation, when Sl^ = Q.. This can be readily understood if we 
consider that the wave energy per unit surface £„ (Eq.l|16ll) can be 
decomposed into two separate terms, as follows; 



Lm^ 161. 

+ -^D(np-n) - 



(21) 



The second term, equal to the angular momentum per unit surface 
times the rotation speed SI, is a local energy transport term (cf. 
Eo.dl lb) and can therefore be represented using the ff-formalism. 
The first term however, equal to the same angular momentum term 
times Qp - Q, is a non local term. In fact the energy flux associ- 
at ed with this non-local transpo rt term is is precisely that identified 
bv lBalbus & Papaloizoijj 19991) as an "anomalous flux", preventing 
self-gravitating discs from acting as pure a-discs. 

We see that far from co-rotation, where Q. t Clp, this term be- 
comes significant, and thus non-local (global) transport becomes 
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important within the disc. For trailing waves launched at co- 
rotation, the direction of energy (and angular momentum) trans- 
port is outwards throughout the disc (i.e. inward travelling waves 
with negative energy density inside co-rotation and outward trav- 
elling waves with positive energy density outside co-rotation), and 
thus the dissipation of such waves effects a net outward transport 
of energy (and angular momentum). If such a wave dissipates at 
large radius (where flp » Q.) then the ratio in which energy and 
angular momentum are added to the disc (Eq.l|20t) is significantly 
greater than the equivalent ratio in the viscous case (Eq.dllb). Con- 
sequently, under such conditions the energy dissipated at large radii 
in a steady state disc with wave transport can significantly exceed 
that dissipated in an equivalent viscous disc ~ with this extra en- 
ergy being extracted by t he wave from deep in the potential and 
transp orted to large radii dLodato & BertirjboOlUBertin & Lodatd 
120011) . However, if waves instead dissipate close to co-rotation, the 
wave transport is dominated by the local term in Eg. (121b: since 
energy and angular momentum transport are exchanged with the 
disc in roughly the same ratio as for a viscous process, then in this 
regime the a-formalism is a good approximation to the actual trans- 
port properties of the disc. 

From Eg. (12 It it is therefore possible to quantify a non-local 
transport fraction ^, from the ratio of the two terms on the RHS, 
such that 



(22) 



Thus in order to assess the importance of non-local effects, we need 
to know the relationship between the angular frequency Q and the 
pattern speed Dp. In section [53] we use the dispersion relation along 
with information extracted through Fourier analysis in order to es- 
timate the pattern speed, and hence to evaluate ^ directly. 



3 SIMULATING THE DISC THERMODYNAMICS 

Realistically simulating the thermodynamics of accretion discs is 
a complex undertaking and as such has received much attention, 
from t he opacity-based treatment employed by Johnson & Gammia 
( |2003|) thr ough t o the various conv e ctive a nd radiative t ransfe r 
models of 'Boss' ('2004'), 'Bolev et al.' ('200?'), 'Maver et al.' llQOf), 
[Stamatellos & Whitworth (2008a) and Stamatellos & Whitworth 
( l2008br) . the latter two of which also account for heating from the 
central star. 

However, in this paper, one of our aims is to investigate the re- 
lationship between the properties of the density perturbations and 
the rate at which the disc cools. This purpose is served most read- 
ily by imposing a known cooling rate against which the heating 
rate and subsequent disc structure may be easily correlated. It is 
therefore not necessary to consider the exact physics of the cool- 
ing regimes found in astrophysical discs, and hence we can use a 
cooling law for the heat loss rate per unit mass QT , such that 



(23) 



where u is the internal energy per unit mass and where the details 
of the cooling function (and possibly of additional external heating) 
are absorbed into the simple parameter f(.ooi- As long as such a char- 
acteristic timescale can be defined, it is therefore possible to use 
this formalism to represent a wide range of cooling mechanisms. 
Within this paper, we use a fixed ratio between the local dynami- 
cal and cooling timescales, such that ji = iltcoo\ is constant. This 
form of cooling has been used extensively in simulations of discs 



in var i ous contexts, for example Gammid j200lh . iLodato & Ricd 
( l2005l) . lHobbs & Navakshii] jlOm . and has proven useful in elu- 
cidating the properties of the gravitational instability in controlled 
numerical experiments. 

In terms of the heating from the gravitational instability, we 
noted in the previous section that density waves extract energy from 
the disc. In addition to compression heating, in the case where the 
pattern speed differs from the rotation speed by more than the local 
sound speed, these waves will steepen into shocks, liberating fur- 
ther heat into the disc. We expect the rate at which energy is added 
to the disc to scale with the energy of the wave and the local dy- 
namical timescale, and we can then express the heating rate per unit 
mass due to the instability as 



(24) 



where we define the radial phase and Doppler shifted phase Mach 
numbers to be Al = |vp|/Cs and M = |Vp|/cs respectively. In Eq. 
( I24t we have also introduced a dimensionless proportionality fac- 
tor e, hereinafter referred to as the heating factor If the relationship 
between the pattern speed and the angular speed is self-similar (i.e., 
it does not vary across the disc), we expect e to be simply a con- 
stant, independent of radius. 

Once the gravitational instability has been instigated and has 
subsequently saturated, we may consider the disc to be in dynamic 
thermal equilibrium, such that the energy released through wave- 
driven shock and compression heating is balanced by the imposed 
cooling, i.e. + = 0. Recalling that u = cl/y(y - 1), we can 
equate equations i24l and i23l and thereby determine the following 
relationship between the amplitude of the density perturbations and 
strength of the cooling, as measured by the /? parameter; 



e/3 7(7- I) \mM 



(25) 



In this paper we shall therefore use global numerical simula- 
tions to test the above energy balance, and to investigate the relative 
magnitude of the local and and non-local transport terms. 



4 NUMERICAL SET-UP 
4.1 The SPH code 

All of the simulations presented hereafter were performed using a 
3D smoothed particle hydrodynamics (SPH) code, a Lagrangian 
hydrodynamics code capable of modelling self-gravity (see for 
example, Benz 1990; Monaghan 1992). All particles evolve ac- 
cording to individual ti me-steps governe d by the Courant condi- 
tion, a force con dition ( lMonagharj|l992h and an integrator limit 
dBate etaljl993) . 

We note here that with SPH, the integral of a physical quantity 
A over a given volume V is estimated by the sum over the individual 
particle values of this quantity, as below; 



X 



AdV.y"-^A., 

VP. 



(26) 



where ra, is the particle mass, and ; loops over all the particles 
within the volume V. In a similar manner, we note that a volume- 
averaged value for A, which we shall call A, can therefore be esti- 
mated via 



(27) 
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when, as in all our simulations, all particles have the same mass. 

We have modelled our disc systems as a single point mass 
(onto which gas particles may accrete if they enter within a 
given sink radiu s, and satisfy certain boundness conditions - see 
iBateetalJfTggl) orbited by 500,000 SPH gas particles; a set up 
common to many other SPH simulations of such systems, (e.g. 
iLodato & Ricel2004l200^ : lRice et alj2003l : IClarke et alj2007h but 
with increaesd resolution. The central object is free to move under 
the gravitational influence of the disc. In order to ensure the sim- 
ulations were properly converged, resolution checks were under- 
taken with discs consisting of both 250,000 and 1,000,000 particles 
- these are discussed briefly in Appendix A. 

As described in section[3]we use a simple cooling model, im- 
plemented in the following manner 

^ = (28) 

where the m; and t^aou are the internal energy per unit mass and the 
cooling time associated with each particle respectively. Again as 
above the functional form of the cooling time is kept simple, such 
that Ojfcooi.i = /3> where Q, is the angular velocity of each particle, 
and where p is held constant throughout any particular simulation. 
All simulations have been run modelling the particles as a perfect 
gas, with the ratio of specific heats y = 5/3, heat addition being 
allowed for via PAY work and shock heating and with the cooling 
implemented as specified above. Artificial viscosity has been in- 
cluded through the standard SPH formalism, with q-sph = 0. 1 and 
ySsPH = 0.2. Note that these values are smaller than those commonly 
used in SPH simulations; we use these valu es to limit the transpor t 
induced by artificial viscosity. As shown in iLodato & Ric3 ( l2004h . 
with this choice of parameters the transport of energy and angu- 
lar momentum due to artificial viscosity is a factor of 10 smaller 
than that due to gravitational perturbations, while we are still able 
to resolve the weak shocks occurring in our simulations. 

By using the cooling prescription outlined above, the rate at 
which the disc cools is governed by the dimensionless parame- 
ter p and the cooling is thus implemented scale free. The govern- 
ing equations of the entire simulation can likewise be recast in di- 
mensionless form. In common with the previous SPH simulations 
mentioned above, we define the unit mass to be that of the cen- 
tral object - the total disc mass and individual particle masses are 
therefore expressed as fractions of the central object mass. We can 
self-consistently define an arbitrary unit (cylindrical) radius Rq, and 
thus, with G = 1, the unit time is the dynamical time f(jy„ = fl"' at 
radius R = \. 
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Table 1. Details of numerical simulations. Note that the duration is quoted 
in terms of the thermal time, equivalent to the cooling time at the outer 
radius « 125/J code units. The j8 = 4 case fragmented, and therefore did not 
run for as long as the other cases. 
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Figure 1. Profiles of Q against radius for different values of the cooling 
parameter fj (top) and mass ratio q (bottom) plotted at the times quoted in 
TableO 



4.2 Initial conditions 

All our simulations model a central point object of unit mass M, = 
1, surrounded by a gaseous disc of mass Mdisc. Although the bulk of 
the simulations have been conducted with a disc to central object 
mass ratio q = M./Mjisc of 0.1, simulations have been run with 
q = 0.05, 0.075 and q = 0.125 to investigate the effects of the mass 
ratio on the non-local energy transport fraction ^. 

All simulations have used an initial mass surface density pro- 
file E ~ which implies that in the marginally stable state 
where Q ~ 1, the disc temperature profile should be approxi- 
mately flat. Since the surface density evolves on the viscous time 
fvisc » 'dyn = ' this profile remains roughly unchanged through- 
out the simulations. The initial temperature profile is ~ iJ"''- and 
is such that the minimum value of the Toomre parameter = 2 
occurs at the outer edge of the disc. In this manner the disc is ini- 



tially gravitationally stable throughout. Note that the disc is not ini- 
tially in thermal equilibrium - heat is not input to the disc until 
gravitational instabilities are initiated. 

Radially the disc extends from R^^ = 0.25 to R„^i = 25.0, as 
measured in the code units described above. The disc is initially 
in approximate hydrostatic equilibrium in a Gaussian distribution 
of particles with scale height H. T he azimuthal v elocities take into 
account both a pressure correction l lLodatoll2007h and the enclosed 
disc mass. In both cases, any variation from dynamical equilibrium 
is washed out on the dynamical timescale. 

Given the dimensions above, one outer dynamical timescale of 
the disc corresponds to 125 time units. To ensure that thermal equi- 
librium is reached and that the gravitational instability is saturated, 
all simulations are followed for at least 10 outer cooling times. To 
this end we shall refer to the thermal time fiherm for each simulation 
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Figure 3. Surface density structures for discs where tlie mass ratio q = with /3 = 5 (left) and /? = 10 (right). The logarithmic scales show surface density 
contours from 10"' to 10"- in code units. Note that the direction of rotation is anticlockwise, and that the plots are given at the times quoted in Table[T] 
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Figure 2. Disc scale height over radius (HjR) plotted as a function of radius 
for varying /3 (top) and q (bottom). Note that in all cases HjR x ^/4 as 
expected. 



as the cooling time evaluated at the initial outer edge of disc, taken 
to be at = 25 - thus fiherm = ?cooi(25) = /3fi"'(25). 



4.3 Simulations run 

In all a total of nine distinct simulations were run for various values 
of the cooling parameter yS and the disc to central object mass ra- 
tio q, as detailed in Table [T] Although previous investigations with 
£ ~ have found tha t the fragmentation boundary is at Ps^^, » 6, 
jRice etalj|200ll2003h . we find that in the case where 1, ~ R ' 
the fragmentation boundary is slightly diflterent, with 4 < ySfma < 5. 



The simulation where /3 = 4 therefore contains a fragment, and 
is included primarily for completeness. All results henceforth are 
given at the time quoted in the final column of Table[T]unless oth- 
erwise stated. The raw data are time-averaged over 500 unit times 
about these values to enhance the signal-to-noise ratio and to give 
the approximate steady-state values. 



5 SIMULATION RESULTS 

Common to all the simulations run is an initial phase in which 
the discs cool rapidly until the value of Q becomes approximately 
unity, at which point the gravitational instability is initiated and 
heat is liberated to balance the cooling. This stage is complete after 
approximately one thermal time, and from then on the discs set- 
tle into a quasi-steady state characterised by the presence of spiral 
arms throughout almost their entire radial range, with Q ~ I. The 
quasi-static Q profiles to which the discs converge are shown in 
Fig. [T] with the cooling parameter j8 varying in the top panel, and 
the disc to central object mass ratio q varying in the bottom panel. 
Note that the data are plotted at the times given in Table[T] Through- 
out all the simulations it can be seen that the discs self-regulate to 
the marginally stable Q ~ I condition over a large range of radii. 

Once the disc has reached a quasi-steady state, the disc as- 
pect ratio H/R also stabilises to the value predicted by the self- 
regulation condition 2 ~ 1 • 



H 

"r 



nI.(R)R^ 



(29) 



and is shown as a function of radius in Fig.|2]for different values of 
P (top panel) and q (bottom panel). 



5.1 Saturation amplitude of the instability 

Once the gravitational instability has been initiated, for the case 
where p < /?frag (i.e. where p = 4) the amplitude of the perturba- 
tions required to balance the cooling rises to the point where non- 
linear effects dominate, leading to the fragmentation of the disc 
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Figure 4. Variation of the relative mass surface density perturbation am- 
plitude i5S/£ with radius for various values of the cooling parameter/?. All 
data plotted at the times shown in Table[T] 



Figure 5. Variation of the radially and azimuthally averaged relative sur- 
face density perturbation amplitude i5S/£ with the inverse cooling parame- 
ter 1 Iji. The radial average is calculated over the range 5 < R < 24. 



into bound objects. In the case where fi > /Jf^g however, the am- 
pMtude increases on the dynamical timescale until the disc reaches 
dynamic thermal equihbrium, at which point the amplitude of the 
surface density fluctuations becomes constant, and the heating they 
provide balances the imposed cooling. This is observed in all our 
simulations where p > 5. 

From the simulations we can now test the prediction for the 
saturation amplitude provided by Eq.l|25ll numerically. Fig.|3]shows 
images of the surface density of the disc for the two cases yS = 5 and 
jS = 10, respectively, where in both cases the mass ratio is q = 0.1. 
It can be seen that, while the overall disc structure remains essen- 
tially constant (as confirmed by a more detailed Fourier analysis, 
see below), the spiral wave amplitude as characterised by the sur- 
face density contrasts appears to decrease with increasing /J. Noting 
that the direction of rotation of the discs is anticlockwise, we find 
that throughout our simulations the waves excited are all trailing 
waves - they all point in opposition to the direction of rotation. 

While the SPH code allows us to conduct a global 3D simu- 
lation of discs with relative ease, it does not permit the direct cal- 
culation of an intrinsically 2D quantity, such as the surface density 
perturbation amplitude Therefore, in order to calculate this 

quantity, we overlay a cylindrical grid on the disc, such that each 
cell contains approximately A^ndgh particles, where N„agh ~ 50 is 
the average number of neighbours within a smoothing kernel for 
our simulations. For each annulus of cells we can calculate an av- 
erage surface density S, and by comparing this to the value calcu- 
lated for each cell within this annulus we can evaluate an annulus 
averaged RMS value for the perturbation amplitude (5S/X. This is 
shown as a function of radius R and the cooling parameter /J in Fig. 

s 

From Fig.|4l it is clear that there is an increasing trend in <5E/E 
with decreasing /? and that furthermore, away from the disc bound- 
aries the saturation amplitude is approximately constant with ra- 
dius. The low values for the perturbation amplitude at small radii 
(R < 5) are probably due to the increased number of particles per 
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Figure 6. Variation of the relative mass surface density perturbation ampli- 
tude i5Z/£ with radius for various values of the disc to central object mass 
ratio q. 

grid cell smoothing out the underlying variation. We can however 
characterise the strength of the perturbation by simply averaging 
6'L/f. over the self-regulated portion of the disc, that we define as 
5 < R < 25 (cf. Fig[TJ. Fig. |5] shows the relation between the 
azimuthally and radially averaged amplitude, which we denote as 
(5S/S>, and the cooling parameter y6. Each point represents a sin- 
gle simulation, while the curve shows our best fit to the data using 
the inverse square root dependence predicted by Eq.(l25t. From the 
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Figure 7. Azimuthal mode amplitudes excited at various radii wliere /? = 4 (top left), /? = 6 (top right), /? = 8 (bottom left) and /3 = 10 (bottom right). 



simulations we therefore obtain the following empirical relation- 
ship, for the case where ^ = 0. 1 



(30) 



In a similar manner we can also calculate the variation in (5S/Z 
with q, which is shown in Fig. [6] We see that the strength of the 
perturbation tends to increase with the mass ratio q, although it is 
clear that this dependence on q is rather less than linear. 



5.2 Fourier analysis: azimuthal structure 

From the simulations we have found empirically that the perturba- 
tion strength (si,/!,^ follows a j8"''^ relationship, as predicted by 
Eq.l|25t. However, this equation also shows a clear dependence on 
the wave modes excited within the disc through the action of the 
gravitational instability. To elucidate this relationship further we 
have therefore conducted a Fourier analysis of the wave modes in 
the disc, a full description of which may be found in Appendix B. 
In this section we therefore consider the effects of both the cooling 
(via /?) and the disc to central object mass ratio q on the excitation 
of the azimuthal m wavenumbers - the next section will describe 
the excitation of the radial wavenumbers. 

We find that in general, whatever the imposed cooling regime, 
for a given mass ratio ^ = 0. 1 the distribution of the azimuthal 
wavenumbers determined by the gravitational instability remains 
approximately constant, with the dominant mode at around m x 5. 
The mode distributions at five radii throughout the disc for the cases 
where j8 = 4, 6, 8 and 10 are shown in Fig.|7] It is clear that the spec- 
tral distribution of the modes shows little variation with radius and 
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Figure 8. Variation of the average azimuthal wavenumber excited as a func- 
tion of wavenumber for/3 = 4-10 where q = OA (top) and as q varies with 
/? = 5 (bottom). 



with the imposed cooling, except that the amplitude of the modes 
decreases as p increases, as expected in view of the decreased per- 
turbation amplitudes seen in Figs [3l i] and ID Additionally, Fig. [8] 
(top panel) shows the variation of the average wavenumber against 
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Figure 9. Surface density structures for discs witli q = 0.05 (left) and ij' = 0. 125 (right), with fS = 5. The logarithmic scale shows mass surface density contours 
from 10"^ to 10"^ in code units. These therefore form adrrect comparison with Fig. [3] where q = 0.1, /} = 5. Once again the direction of rotation of the discs 
is anticlockwise. 





Figure 10. Radial mode amplitudes excited at various radii where yS = 4 (top left), /} = 6 (top right), fS = 8 (bottom left) and /3 = 10 (bottom right). 



radius for all the values of /? that have been tested, and we note that 
ahhough some small variation is seen, it is uncorrelated with the 
imposed cooling. 

The bottom panel of Fig. [8] shows the variation in the aver- 
age azimuthal modes excited as the disc to central object mass ra- 
tio q varies, while the cooling is held constant with /? = 5. It can 



be seen from the plot that variation of this parameter does have a 
marked effect on the power spectrum of the waves - the average 
mode number varies inversely with the mass ratio, from m„ « 15 
where q = 0.05 to /Wav ~ 10 where q = 0.125. This variation is 
also clearly seen in Fig.[9l where a large number of flocculent arms 
are present in the disc with q = 0.05, and fewer, rather more well- 
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Figure 11. Mode amplitudes plotted against the product of the peak radial modenumber and the disc scale height H excited for various radii where /? = 4 (top 
left), /? = 6 (top right), /3 = 8 (bottom left) and /3 = 10 (bottom right). 



defined spiral arms a ppear in the disc wher e g = 0. 125 (cf. the simi- 
lar result obtained in iLodato & Ricel2004 . By comparison, the left 
panel of Fig.[3]shows a disc with y6 = 5 and ^ = 0. 1, and the pattern 
of spiral arms present is intermediate to those shown in Fig.|9] 

5.3 Fourier analysis: radial structure 

We now consider the radial wavenumbers k of the waves excited by 
the gravitational instability. In contrast to the azimuthal modes, it 
is clear from Figs. [3] and |9] that there is significant variation in the 
radial wavenumber k with radius, and Fig. |9] suggests that there is 
an additional variation with the disc to central object mass ratio. 

In Fig. [To] we show the variation in the power spectrum of dif- 
ferent radial wavenumbers for the cases where p = 4, 6, 8 and 10 
and ^ = 0. 1, at the same radii as the azimuthal wavenumbers shown 
in Fig. [7] As with the azimuthal modes, we see little overall change 
in the spectral distribution of the modes with varying yS excepting 
that the amplitudes of the modes decrease as the cooling weakens. 
Conversely however, these plots show a significant variation with 
radius, in that the peak wavenumber decreases with increasing ra- 
dius, and thus the dominant wavelength similarly increases with 
radius. 

Figure [TT] on the other hand shows the power spectrum as a 
function of kH, where we normalize the wavenumber to the ex- 
pected most unstable one, (see Section 2.1). We thus con- 
firm the expectations from the linear WKB approach, as these plots 
show a clear peak at kH « 1. This therefore suggests that through- 
out the disc the excited waves are close to co-rotation. Fig.ll2l(top 
panel) shows the average radial wavenumber as a function of radius 
for all the values of /? considered, again confirming the trends al- 



ready discussed and also further showing that, excepting the varia- 
tion in amplitude discussed above, the structure excited by the sim- 
ulation is essentially independent of the cooling imposed. It also 
shows that the simulation to simulation scatter is very small. 

The bottom panel of Fig. [T2] shows that, as with the az- 
imuthal wavenumbers, there is clear variation in the average ra- 
dial wavenumber k^y with the disc to central object mass ratio for 
a given p (in this case p = 5); increasing the mass ratio decreases 
the average wavenumber in approximate inverse proportion. The 
dashed line in the bottom panel of Fig. [T2] plots a simple R^^'^ 
curve, indicating that the average wavenumber follows a power-law 
distribution with radius, such that k.^, ~ R^^'^, which remains con- 
stant with varying mass ratio. This is easily understood by noting 
that since the sound speed c, is approximately constant by construc- 
tion, Eq.lO indicates that k ~ 1, ~ R"^^^. 

5.4 Mach number of the spiral modes 

Returning briefly to the dispersion relation given in Eqs.l[T) and 
l|2j, we note that this is only strictly valid for infinitesimally thin 
discs. As our simulations are fully three dimensional, we require 
a correction to the self-g ravity term t o account for this, and this is 
shown in Eq.l|3l) below ( lBertinll2000l) . 

m^(fip - n)^ = cjfc2 - ^^^^ + n^, (31) 

" 1 + \k\H 

where we have also used the fact that our discs are approximately 
Keplerian, and thus k x Q.. The reduction factor of 1/(1 -I- \k\H) 
arises from the vertical dilution o f the gravitat i onal potential due to 
the finite thickness H of the disc l lBertinl2000l : lBinnev & Tremaind 
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Figure 12. Variation of the average radial wavenumber as a function of 
radius for /3 = 4 - 10 and ^ = 0. 1 (top) and as q varies with /? = 5 (bottom). 



Figure 13. Wave phase Mach number M (thick lines) and the Doppler- 
shifted phase Mach number M (thin lines) as a function of radius for various 
values of /? with ^ = 0.1 (top) and as q varies with yS = 5 (bottom). 



|2008| ; I Vandei^oortll 1 970h . Using this finite-thickness dispersion re- 
lation and our averaged values for k and m it is possible to calculate 
a Doppler-shifted angular speed IDp - n|, noting that the sign of 
flp - n cannot be determined from Eg.Olt. Since the average ra- 
dial wave-number is always very close to the most unstable one, 
and since the disc is almost exactly marginally stable, the resul- 
tant average pattern speed turns out to be always very close to co- 
rotation (as can be seen in Fig. I16t. We quantify the deviation of 
the pattern speed from co-rotation later in section lSTS] 

We can further calculate the radial phase and Doppler-shifted 
phase Mach numbers, and these are shown in Fig. [O] The upper 
panel shows the wave phase Mach number M (thick lines) and the 
Doppler shifted phase Mach number M (thin lines) as functions of 
radius for various values of [} with q = 0.1. Similarly, the lower 
panel of Fig. [T3] shows the variation of these Mach numbers with 
the mass ratio q for - 5. We see immediately that both quantities 
are independent of the cooling rate as measured by fi with very little 
scatter Moreover, the Doppler-shifted phase Mach number is very 
close to unity. In a similar manner this quantity remains unchanged 
with variations in the mass ratio, although the phase Mach number 
decreases with increasing q. 

The above results essentially imply that the wave structure is 
determined by the requirement that the normal component of the 
flow into the shock is almost exactly sonic - a natural criterion for 
a quasi-steady system due to the dissipative nature of shocks. For 
waves with winding angle i, and radial Doppler-shifted phase speed 
Vp, a sonic normal component of velocity into the shock implies 



: Cs, leading to 



1 

cos i 



(32) 



Hence, in the limit of tightly wound waves where cos / « 1, we 
should expect that M. ~ 1, as indeed we find in Figure [T3] For 
completeness. Fig. [14] shows the winding angle ; as a function of 
radius for varying yS (top) and mass ratio q, (bottom), using the 
definition tan; = mIkR. In all cases, ; < 15°, so the waves are 
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Figure 14. Wave winding angle (' as a function of radius plotted against 
varying cooling (top) and mass ratio (bottom). 



reasonably tightly wound throughout. Again there is no significant 
variation with cooling, but the structure becomes more open as the 
mass ratio increases, as expected from Figs.[3]and|9] 

We can also use Eq.l|25|l to estimate the amount of energy 
dissipated by the weak spiral shocks per dynamical time, as char- 
acterised by 6, which is shown in Fig. [15] We have seen through 
the constancy of the Doppler-shifted phase Mach number that the 
shock structure that forms in the disc is indeed self-similar, and 
thus the heating factor 6 is also largely independent of the applied 
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Figure 15. The heating factor e as a function of radius for various values of 
/3 with ^ = 0. 1 (top) and as q varies with /? = 5 (bottom). 



Figure 16. Non-local transport fraction f as a function of radius for various 
values of fi with ^ = 0. 1 (top) and as q varies with /3 = 5 (bottom). 



cooling, the mass ratio and the radial position. Note that the larger 
values for e generated at low radii {R < 5) are probably due to the 
inaccuracies in calculating i5E/E in this region rather than a break- 
down in self-similarity. 

5.5 On the locality of transport induced by self-gravity 

In the previous subsection we noted that the (spectrally averaged) 
pattern speed of the waves f2p is always very close to the angular 
velocity of the flow f2, thereby indicating that the waves are close 
to the co-rotation resonance as suggested earlier in the results of 
the radial mode decomposition. We can estimate more quantita- 
tively how close to co-rotation the spiral waves lie by calculating 
the quantity as given in Ea.ll22t. This is shown in Fig. [16] as a 
function of radius for all the values of j8 and mass ratio simulated. 

For the ^ = 0.1 case, we thus see that varying the cooling 
has no significant efi'ect on the transport properties of the disc, with 
^ » 0. 1 throughout the radial range, albeit with some scatter. This 
means that in this configuration the disc is dominated by local trans- 
port processes, and is as such reasonably well des cribed by the vis- 
cous a prescription of IShakura & SunvaevI ( Il973h - global effects, 
although not negligible, are smaller than local efi^ects by an order 
of magnitude. 

By varying the mass ratio, we see that the strength of non- 
local effects increases with q, rising to f ~ 15% for the case where 
q = 0.125. This confirms the results of lLodato & Ricd | |2004|) . who 
found similarly that non-local effects (characterised by strong tran- 
sient structures in the disc) become increasingly important as the 
disc mass ratio rises, although for the parameter range considered 
here the disc remains dominated by local effects. This non-local be- 
haviour can be elucidated further by noting that from the definitions 
of f and Vp (Eqs.llH and GDl) we have 



m m cos i 



(34) 



\ m 



which, for kH x 1 reduces to 



(33) 



The non-locality of the transport is therefore directly linked to the 
openness of the structure that is induced in the disc through self- 
gravity. Since, as we have noted earlier, for larger disc masses the 
spiral structure tends to become more open and more dominated by 
low m modes, we then see that more massive discs tend to become 
more subject to non-local effects. 



6 DISCUSSION AND CONCLUSIONS 

In this paper we have undertaken 3D global numerical simulations 
of gaseous, non-magnetised discs, evolving under the influence of 
a massive central object and their own self-gravity. We have mod- 
elled the gas as an ideal gas with y = 5/3, together with a simple 
cooling prescription based on a local cooling timescale. We have 
used these simulations to investigate the structure that forms once 
the discs have settled into a quasi-steady marginally stable state as 
a function of both the imposed cooling and the disc to central object 
mass ratio. 

We have found that the amplitude of spiral arms induced in 
self-gravitating discs, as characterised by the RMS surface density 
perturbations, can be described straight-forwardly through the em- 
pirical relationship 



LO 

V8' 



(35) 



(where fi is the ratio between the local cooling and dynamical 
timescales), with only a weak dependence on the disc to central 
object mass ratio. This is in fact closely linked to the result that 
the Doppler-shited Mach number is very close to unity - by con- 
sidering the entropy change AS across an adiabatic shock where 
the Mach number Ma 1, it can be shown that A5 ~ (M^ - 1)^- 
Thermal equilibrium in these discs is established between cool- 
ing, at a rate inversely proportional to /?, and the irreversible con- 
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version of mechanical energy into heat, at a rate proportional to 
the entropy jump A5 at the shock front - we therefore find that 
[} ~ (M- - . Standard shock relations show that the density per- 
turbation 6plp ~ {M- - 1), and hence simply from considering the 
properties of weak adiabatic shocks we can arrive at the relation- 
ship 5p/p ~ /S^''^ 

Additionally we find that the heating factor e - that fraction of 
the available wave energy that is liberated as heat back into the disc 
gas - remains essentially invariant at x 20% with both the imposed 
cooling regime and the mass ratio of the disc to the central object. 

As expected, our simulations show that the dominant radial 
wavenumber is approximately equal to the reciprocal of the local 
scale height of the disc throughout the radial range, k x nCL/cl. 
We therefore find that the radial spacing of the arms is dependent 
only on the surface density and temperature profiles of the disc. 
Likewise although further work is required to understand the rela- 
tionship fully, the azimuthal disc structure is dependent on the disc 
to central object mass ratio, with more massive discs being charac- 
terised by more open structures than their lower mass counterparts 
for a given central object mass. 

Our numerical r esults bear out the t heoretical analysis of 
iBalbus & Papaloizoij ( Il999h and iGammid 1^01), who suggest 
that discs in the Q ~ I marginally stable state may be mod- 
elled as predominantly lo cal. Simulat i ons of self-gravitating discs 
with radiative transfer bv lBolev et alj j2006h also found that close 
to co-rotation, angular momentum transport was well modelled 
by a local a-prescription ev en when global modes were present. 
iBalbus & Papaloizoul [m^ further predicted that non-local trans- 
port from an "anomalous flux" proportional to Q - flp would be- 
come significant far from co-rotation, a result we have derived an- 
alytically using the WKB approximation for tightly wound waves. 
We have then used the WKB dispersion relation along with em- 
pirically determined information on the dominant wavenumbers to 
make an estimation of l^p - Q.\. We find that, at least for low mass 
discs, this is a small fraction of Q (less than 15% for discs with 
q < 0.125, regardless of the efficacy of the cooling). Our results on 
the magnitude of the non-local transport fraction ^ = \Q.p - n\/n 
can furthermore be readily understood in terms of the empiri- 
cal constancy of the Doppler-shifted radial phase Mach number, 
M. We conclude that the importance of such non-local effects in 
gaseous self-gravitating discs is set by the self-adjustment of the 
pattern speed to ensure that the normal flow speed into the arms is 
sonic. We have then demonstrated that this condition implies that 
^ » m"' sec i, where / is the opening angle of the spiral structure. 
Since the structure within the disc becomes more open as the disc to 
central object mass ratio increases, this implies that the importance 
of non-local transport also scales with q. 

We note also that in coUisionless sytems such as stellar discs, 
this self-regulation process for the pattern speed breaks down as 
shocks cannot form. Hence it is possible to excite global modes in 
such discs, and thus non-local transport of energy and angular mo- 
mentum may be more significant dynamically. The results that we 
present here are therefore restricted to the case of predominantly 
coUisional, gaseous discs. Our results provide a theo retical under- 
pinning for the results of jLodato & Ric3l2004[2005il on how the 
importance of global transport depends on the disc to central object 
mass ratio in gaseous discs. In particular, we note that in cases (like 
those described here) where the disc mass is a small fraction of 
the central object mass (as could be the case for relatively evolved 
protostellar discs) the effects of self-gravity are expected to be well 
described as a pseudo- viscous process. 

One of the most important applications of our study is that 



we can relate the amplitude of spiral modes in gaseous discs to the 
cooling regime. With ALMA coming online in the relatively near 
future, promising milli-arcsecond resolution in the millimetre/sub- 
mm range, it is possible that such observations of spiral structure in 
proto-planetary discs may become technically feasible. 
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Figure 17. Variation of the radial RMS surface density perturbation ampli- 
tude as a function of simulation resolution. 



APPENDIX A: RESOLUTION AND CONVERGENCE 
TESTS 

In this appendix we shall briefly outline the tests that were under- 
taken to ensure the convergence of our results. 

Three simulations were run, all with the cooling parameter 
j8 = 6 and mass ratio q = Q.l, using discs of 250,000, 500,000 and 
1,000,000 particles. These were otherwise identical to the simula- 
tions that were used for this paper, as described in full in section 
O The three values that are of most significance to our results are 
the RMS surface density perturbation amplitude <5Z/S and the av- 
erage radial and azimuthal wavenumbers k^^ and niav respectively. 
Fig.|17|shows how (5X/£ varies with resolution, and although there 
is considerable scatter it is clear that there is no systematic varia- 
tion with resolution. A similar result (with even less scatter) is also 
obtained when one conducts a Fourier analysis of the simulations - 
there is no systematic variation with resolution. We may therefore 
conclude that our simulations are converged, and that the resolution 
when using 500,000 particles is satisfactory for our purposes. 



APPENDIX B: FOURIER DECOMPOSITION METHODS 

In this appendix we detail how the Fourier mode analysis was con- 
ducted, using the SPH particle positions as the input values. For 
simplicity, we begin by discussing how the radial k mode ampli- 
tudes were computed, as this had practical implications on how the 
azimuthal m mode analysis. 



Radial mode analysis 

To calculate the radial Fourier mode amplitudes within a disc 
presents certain problems, since even a cursory glance at Fig.|3]re- 
veals that the radial wavenumber k varies significantly with radius. 
Also, unlike the azimuthal wavenumbers, the disc is neither uni- 
form nor periodic in radius, and therefore the underlying Fourier 
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Table 2. Details of the Fourier analyses for the various disc to central object 
mass ratios analysed. 



distribution corresponding to the disc surface density profile also 
has to be taken into account. The following method addresses both 
of these problems while keeping the signal to noise ratio as high as 
possible. 

The disc to be analysed is divided into numerous overlapping 
annuli of width hR, which varies with the central radius of the an- 
nulus, and into a number of sectors of fixed angular width Ai^. The 
Ai? values are chosen such that each annulus is of sufficient ra- 
dial extent to resolve the greatest radial wavelength present at that 
radius, likewise each sector must be narrow enough to ensure the 
wave crests are distinct and not smeared out across a wide range in 
R. In this manner, the smaller the winding angle = tan"' \mlkR\ 
of the waves the wider the sectors can be for a given resolution. 

Since the radial wavenumber profile depends on the disc to 
central object mass ratio, the radial extent /S.R of the annuli varied 
likewise, in order to capture all the relevant modes. The values used 
in our analyses are summarised in Table[2] Note that the widths of 
the annuli increase linearly across the disc, from the initial to the 
final widths quoted. 

To calculate the underlying Fourier distribution due to the un- 
perturbed surface density profile, the Fourier transform was taken 
over the whole of each annulus. This thereby smears out all the 
waves and takes the average distribution, and is evaluated accord- 
ing to the following relation; 



1 

N 

^ ' ann 



(A-1) 



where is the k mode amplitude corresponding to the underlying 
disc distribution, A'ann is the number of particles per annulus, k is 
the radial wavenumber and the Rj are the radii of the individual 
particles. 

The Fourier distribution of the waves overlaid on the disc are 
calculated by taking an equivalent transformation over each sector 
within the annulus, such that 



1 

A'sect 



(A-2) 



where Ai; „ is the k mode amplitude of the waves and disc evaluated 
in the nth sector, and Nsea is the number of particles in that sector. 

Finally the Fourier distribution due solely to the waves in each 
sector is given by the difference between equations iA-2\ and ( lA-lb . 
Since each sector should be statistically similar to the others we 
may then average over all the sectors Nseaoa (which in this case 
does not smear the wave component out, but reduces computational 
noise), to give the average radial Fourier mode amplitudes of the 
waves {Ak); 



<A,>= Yj (^tn-^k)- 



(A-3) 
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Figure 18. Test case for Fourier analysis showing the imposed structure. 
The colour scale shows logarithmic surface density. 
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Azimuthal mode analysis 

For the azimuthal m wavenumbers the analysis is more straightfor- 
ward. The disc was initially divided into annuli of fixed width /S.R 
in such a manner that each of these annuli is narrow enough to en- 
sure the wave crests occupy only a small range in (f>. In contrast to 
the radial modes, these annuli therefore need to become narrower 
with decreasing winding angle to maintain resolution. We found 
AW = 0.2 (in code units) to be sufficient for the purposes of this 
analysis. The azimuthal wavenumber amplitudes A„, within each 
annulus are then computed via 



Figure 19. Results of the Fourier decomposition of test disc shown in Fig 
llSl in terms of the azimuthal wavenumber, m. 
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^ Linn 



-im<Pj 



(A-4) 



where the 0, are the azimuthal angles of the individual particles, 
A^ann the number of particles in each annulus and m the radial 
wavenumber of the wave, corresponding to the number of arms in 
the spiral. 

However, to ensure that we have the azimuthal m-mode am- 
plitudes specified at the same radii as the radial A:-modes, then an 
average value is taken of the ni-mode amplitudes over all annuli 
where the central radius falls within that annulus in which the k- 
modes are determined. 



Analysis Checks 

To ensure that the results of the Fourier analysis are accurate, we 
ran the following test case. A disc with an underlying surface den- 
sity profile E ~ R^^^^ and an analytically superimposed structure 
was created with five spiral arms, such that m = 5 throughout 
the entire disc and with the radial wavelength increasing linearly 
from /!„,;„ = 2 to A^^^ = 7.6. This gives a total of five full wave- 
lengths across the face of the disc, which extends from = 1 to 
R = 25. The surface density of the disc, clearly indicating the im- 
posed structure, is shown in Fig.[T8] The Fourier analysis was then 
conducted using the annulus widths quoted in Table [3] where as 
described above the width of the annuli increased linearly from the 
minimum to the maximum quoted value. 

The results from the azimuthal Fourier decomposition are 
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Table 3. Details of the Fourier analyses for the test case 



shown in Fig. [191 ™d show that the azimuthal wavenumber is re- 
solved extremely well. The fundamental frequency m = 5 is clearly 
dominant, with no other modes except higher harmonics present at 
any significant amplitude. Note that the results for the azimuthal 
modes show no sensitivity to the annuli used for the analysis, and 
are quoted for the first case in Table [3] where the annuli width cor- 
respond exactly to the radial wavelengths. 

The results of the Fourier decomposition for the radial 
wavenumbers are shown in Fig. [20] which shows the actual distri- 
bution of wavenumbers (as calculated directly from the known dis- 
tribution of wavelengths) and the distributions derived from analy- 
ses using the annuli given in Table [3] Note that the analysis using 
annuli that fit the wavelengths exactly correspondingly reproduces 
the exact result. Clearly there is scatter within the results for the 
radial wavenumbers, which arises from the fact that if the actual 
wavelength is not an integer divisor of the annulus over which the 
analysis is being conducted, more than one wavenumber appears to 
be excited. We note however that the scatter is never more than a 
factor of 1.5 above or below the true value, which we deem to be 
accurate enough for the purposes of this analysis. 
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Figure 20. Results of the Fourier decomposition of tlie test disc sliown in 
Fig. llSl sliowing the peak radial wavenumber <:„,ax as a function of radius. 
Various analyses are shown using the annuli given in Table[3] 



Figure 21. Resolution limits for the Fourier analysis in terms of the az- 
imuthal wavenumbers (top) and radial wavenumbers (bottom). The dashed 
line in the bottom plot indicates /cmax = 10. The simulation parameters are 
/?= 10,^ = 0.1. 



Resolution Limits 

Throughout this paper, we have used discs of 500,000 particles 
when undertaking the Fourier analysis. From the fundamental SPH 
resolution limit of the local smoothing length h we can evaluate the 
maximum resolvable azimuthal and radial wavenumbers m^^.^^ and 
^n,ax as a function of radius throughout our discs, such that 

2nR lR\lH\ 
.W = — = 2. (-)(-), (A-5) 

In 27t/H\ 

Since the approximate expected values for radial and azimuthal 
wavenumbers are such that kH ^ 1 and mH/R a 1 respectively, 
EQS.l lA-5t and ( IA-6I 1 show that the accuracy of the Fourier analy- 
sis is closely tied to the vertical resolution of the disc through H/h. 
Using the average smoothing length at each radius, these resolution 
limits are shown in Fig.[2T| We have used data from the simulation 
where = 10, as this gives the most conservative limits of all our 
experiments. 

The vertical resolution of the disc as indicated by H/h is 
shown in Fig. [22] for simulations using 500,000 particles. We find 
that the disc thickness is covered by approximately two smooth- 
ing lengths throughout, and thus is adequately resolved. For the 
Fourier analysis we therefore see that the expected peak wavenum- 
bers are resolved by a factor of approximately 4n throughout the ra- 
dial range. For the radial wavenumbers, we are primarily interested 
in A; < 10, which is well resolved until at least R = 25, again ade- 
quate for the analyses we have undertaken. We conclude therefore 
that throughout the radial ranges of interest, the Fourier analyses 
we have presented are well resolved. 
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Figure 22. Ratio of the disc scale thickness H to the average smoothing 
length h as a function of radius. Again the simulation parameters are /? = 10, 
q = 0.1. 



